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Abstract. We describe and discuss the properties of three numerical methods for solving the diffusion equation 
for the transport of the chemical species and of the angular momentum in stellar interiors. We study through 
numerical experiments both their accuracy and their ability to provide physical solutions. On the basis of new 
tests and analyses applied to the stellar astrophysical context, we show that the most robust method to follow the 
secular evolution is the implicit finite differences method. The importance of correctly estimating the diffusion 
coefficient between mesh points is emphasized and a procedure for estimating the average diffusion coefficient 
between a convective and a radiative zone is described. 
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1. Introduction 

For many physical processes in stellar interiors as heat 
transfer, transport of angular momentum, mixing of chem- 
ical elements through time dependent convection, semi- 
convection or turbulence, a diffusion equation needs to 
be solved. Diffusion is thus involved in numerous inter- 
esting astrophysical problems. For instance, Vauclair and 
Charbonnel (19-98,) discuss the importance of selective 
diffusion in low-metallicity stars and the consequences 
for lithium primordial abundance; Charbonnel l|1995fl has 
shown how rotational diffusion in low mass stars may lead 
to the destruction of ^He and thus to new insights on the 
evolution of the abundance of this cosmological element; 
most of the results on the effects of rotation in stellar in- 
teriors are based on the resolution of a diffusion equation 
for both the transport of the chemical species and of the 
angular momentum (Chaboyer & Zahn lTM^ Zahn 1992 
Talon & Zahn lTMTI Denissenkov et aL lTMHI Heger et al. 
EnnOlMaeder & Mevnet IMTTl Meynet & Maeder EOHI ) . 
Diffusion plays also a key role in our understanding of 
the structure and the cooling of white dwarfs (cf. Kawaler 

For solving numerically the diffusion equation, differ- 
ent methods are available. Quite generally, they can be 
classed into two main categories: the finite differences 
methods and the finite elements methods. Depending on 



how the time discretisation is performed, one distinguishes 
three subclasses in each of these two main categories, 
namely the subclasses of the explicit, implicit and of the 
"Crank-Nicholson" type methods. Thus, at least six differ- 
ent methods are available. Here we are interested in mod- 
elling the long term evolution of stars (secular evolution). 
This immediately prevents us of using explicit methods 
which requires the use of much too short time steps (see 
below). Among the four remaining methods, the "Crank- 
Nicholson" finite elements method was not tested in view 
of the results obtained by the "Crank-Nicholson" finite 
differences method (see Sect. 5.4). Thus we shall focus 
our attention on three of them namely 

— the implicit finite elements method, 

— the "Crank-Nicholson" finite differences method, 

— the implicit finite differences method. 

Descriptions of these methods as well as a general dis- 
cussion of their respective advantages and weaknesses can 
be found for instance in Press et al. I|1992|l for the finite 
differences methods and in Zienkiewicz & Taylor H2()()()|l 
for the finite elements methods. If these considerations can 
help in choosing the best method in general, we think that 
very valuable complementary information can be gained 
by applying the methods to the resolution of realistic as- 
trophysical cases, where the different timescales involved 
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are very specific. Tliis is wliat we propose liere and tliis is 
the main aim of this work. Let us recall that, as said in 
Press et al (|iy921 "... differencing partial differential equa- 
tions is an art as much as a science", and that a numerical 
scheme, which from a theoretical point of view does appear 
as very robust can show some weaknesses when applied to 
a real particular situation. This is precisely why it is nec- 
essary to test realistic astrophysical situations. This allows 
us to make a new and appropriate analysis to support the 
choice of the best method to handle diffusion problems. 

The interested reader will also find in this paper a de- 
tailed description of the discretized equations for resolving 
the diffusion of chemical elements and of the angular mo- 
mentum in stars. We also propose a new way to estimate 
the diffusion coefficient in the interface between a convec- 
tive and a radiative zone. 

In Sect. 2 we briefly recall the physical problem to be 
resolved. The three numerical methods are described in 
details in Sect. 3. The estimate of the diffusion coefficient 
at the interface between a convective and a radiative zone 
is presented in Sect. 4. Sect. 5 is devoted to various tests 
and comparisons. Section 6 summarizes the main results. 

2. The physical problem 

Diffusion is a process by which components of a mixture 
move from one part of a system to another as a result 
of random motion. Let us briefly recall how the diffusion 
velocity is defined (see e.g. Battaner I1996f) . In a multi- 
component fluid, one defines the mean velocity Vq of the 
mixture as 



Vo = 



The mean value of Vj is 



< V, >= — I f^V^dTp =< V, > -Va. 



(1) 



This quantity is called the diffusion velocity of the ith 
component. It can be shown very easily, using the expres- 
sion for that the different diffusion velocities must com- 
pensate one another, i.e. 



mini < Vi >= 0. 



(2) 



There is thus no net mass flux associated to diffusion. In 
that respect, in a star, the changes of chemical compo- 
sition due to diffusion can be treated very similarly as 
those due to the nuclear reactions. Like the nuclear reac- 
tions, diffusion modifies the mean molecular weight (and 
thus the hydrostatic structure of the star), but does not 
induce any net mass flux. 

In the context of the Boltzmann microscopic theory, it 
is possible to deduce from the equations of motion for the 
different components, expressions for the < Vi >'s (see 
Battaner 1996). We shall not repeat these developments 
here. Instead we consider that diffusion is equivalent to a 
displacement of particles whose velocities may be related 
to the diffusive coefficient, D, by the expression (see also 
Sect. 3.2) 



< V, >-- 



D 



where 



where Xi is the mass fraction of element i. The minus 
sign results from the fact that the velocity is directed in 
the opposite direction to the mass fraction gradient. This 
expression automatically satisfies the condition seen above 
that '^irmrii < Vi >= 0, indeed 



<Vi>^ — / fiV.dTp, 
Jp 

is the mean velocity of the ith component, Ui is the par- 
tial number density of particles i, mi their mass, Vi their 
velocity and fi(x,p,t) the probability distribution func- 
tion, which gives the probability to find a particle i at 
the position x, with the momentum p, at the time t; 
Tp is the phase space volume for the particle's momen- 
tum. One can wonder, why in the expression for Vq, the 
weighting factor is the partial mass density and not the 
partial number density. The question could be formulated 
in another way: when is the mean velocity of the mixture 
zero ? When the net flux of the number of particles is 
zero C^i rii < Vi >— 0), or when the net flux of the mass 
is zero (^^ mini < Vi >= 0) ? Obviously when the net 
flux of the mass is zero, otherwise the centre of gravity of 
the system is moving and thus the mean velocity of the 
mixture is not zero. 

The peculiar velocity of a particle is deflned as 

Vi = - Vq. 



since Xi — I and the relation between n; and Xi is 

pX, 

ni = , 

AiWiH 

where p is the density, Ai the atomic weight expressed in 
units of proton mass mu . 

We are interested in resolving the diffusion equation 
in spherically symmetric systems. Let us consider a one- 
dimensional problem, where particles i may diffuse along 
the radial direction r. Thus 



< y^ > -- 



D dXi 

x'~d^ 



(3) 



The variation of the number of particles in the element of 
volume V and of surface S, due only to diffusion, is 



^ J nidv = - J n^ < Vi > -dS. 



(4) 
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From the divergence theorem and the expression of < Vi > 
just seen above one obtains 



dt 



j,2 Qj. 



{r'^n, <Vi>) 




d 



pr^D- 



(z) 

dr 



.(5) 



Replacing rii by its expression as a function of Xi in Eq. (jSJl 
gives 



d_ 

dt 



1 d ^ 2jj'^Xi ^ 



(6) 



Thus Eq. © results from the expressions for the <Vi> 
(Eq. O and from the continuity equation for the n^'s 
(Eq.EI). 

Sometimes, Eq. ((BJ is written as below 



d Id 2r^^C^ 



(7) 



with partial concentration Ci instead of mass fraction, 
without indication on the precise meaning of what con- 
centration means. Partial concentration in mass is equiv- 
alent to the mass fraction, while partial concentration in 
number is equal to rii/n where n — X^i"-*- Replacing Ci 
in Eq. (|7|) , by /n and expressing rii as a function of Xi , 
one obtains 



dt 



dr 



where p, — p/{nmu) is the mean molecular weight of the 
ions. This equation is identical to Eq. only when p 
remains constant as a function time and is constant as a 
function of r. Thus Eq. ||SJ| is equivalent to Eq. © only 
when applied to minor constituents, the abundances of 
which do not affect /i and when p has no gradient. In 
more general cases, Eq. © is not equivalent to Eqs. (01 
and (0. 

This can also be seen in a slightly different way. When, 
m Eq. ((7Il,thec,'s are identified with the number fractions, 
one obtains 



d / ni\ i C / 2 "-i 
dt\n) r'^ dr \ ^ n 



1__9 

J.2 Qj, 



(9) 



where we have introduced a velocity Wi — ^--^ (—)■ 

^ * iii I n or \ n I 

This velocity is not a diffusive velocity in the sense that 
'^iTniniWi — —nDnindp,/ dr is not equal to zero except 
in zones where p, is constant. Eq. jnj expresses the change 
of prLi/n = pniTTiH in an element of volume resulting from 
the transport of the quantity pui/n. But physically, this is 
not prii/n which diffuses but the particles themselves. As 
above, there are two conditions for Eq. Q to be equiv- 
alent to Eqs. (jSJ and 1) the particles which diffuse 
must have a very small abundance, so small that their 
diffusion does not affect p = p/{nmH)', 2) and p does 
not vary with the radius. We conclude thus that when 
the diffusive transport is described by an equation of the 
form given in Eq. (0, the concentrations are equivalent to 



mass fractions. The identification with number fractions 
results in unphysical description of the process except in 
very particular situations. 

Starting from Eq. ^ and summing over all the chem- 
ical species i, one obtains 



dp I d ( .d 



dt 



dr 



dr 



1' 



0. 



This is quite consistent with the fact that the diffusive 
velocities must compensate each other. Thus the density 
can be put outside from the time derivative in Eq. 



dX, 



j,2 Qj, 



pr^D 



dX, 
dr 



(10) 



where is the langrangian mass coordinate. This is the 
equation we shall numerically resolve in this paper. The 
conditions at the centre and the surface are: 



dX, 
dr 



dX, 
dr 



= 0, 



(11) 



where M is the total mass of the star. 

The equation expressing the diffusion of the angular 
momentum is (see e.g. Endal & Sofia IIOTT^ 



(8) d{r^Vt) 



dt 



1 5_ / 4^/917 

J.2 Qj, K Qj, 



(12) 



where $7 is the angular velocity and D' the diffusion co- 
efficient for the angular momentum. In our numerical ex- 
periments, the radii do not change with time^. Eq. H12|l 

then becomes 



pr 



'dt 



rrii 



j,2 Qj, 



= div {pr^D'V{n)) .(13) 



The conditions at the centre and at the surface are 

= 0. (14) 



dn 

dr 



dn 

dr 



=M 



Let us notice that in general the transport equation for 
the angular momentum may contain other terms express- 
ing the advection of angular momentum by meridional 
circulation and the effects of magnetic braking. As these 
last two terms bring no specific problems with respect to 
diffusion, we do not consider them here. However, the ad- 
vection of angular momentum needs a particularly careful 
treatment as well (cf. Meynet & Maeder 2000). 



^ Let us note that in more realistic stellar models, the radii 
of course vary as a function of time, however time steps can be 
chosen sufficiently small for considering them to be constant 
during one time step. The same can be said for the diffusion 
coefficient and other structure variables as for instance the den- 
sity. 
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3. The three tested methods 

3.1. The implicit finite elements method 

The detailed description of the imphcit finite elements 
method used to resolve the diffusion equation for the 
chemical species is presented in Schatzman et al. p981|l 
(see the appendix by Glowinsky & Angrand). We present 
here the procedure for the case of the angular momentum 
diffusion. The basic idea of this method is to decompose 
the unknown function, here fl(rnr,t), as a linear combi- 
nation of well chosen independent functions. This decom- 
position can be performed in many different ways. We 
adopt here the same decomposition as in Schatzman et al. 
H1981|l . Of course, the conclusions concerning the ability 
of this method to provide physical solutions will only refer 
to this particular choice. 

We decompose the star in K shells, with the lan- 
grangian mass coordinate of the i*'* mesh point being rrii 
{rrii is the mass inside the sphere of radius Vi). In the fol- 
lowing all the quantities with an indice i are evaluated at 
the mesh point i. The shells are numbered from 1 at the 
surface to K at the centre. Let us introduce K functions 
bi{mr) defined by 



— II rrir G TO.i,m,_i , 



h{mr) 



\ 







if TUr ^ [rrii+i, rrii 



The function bi is equal to one at = rrii, is equal to 
zero at mi+i and mi_i and varies linearly as a function 
of rrir inbetween. By multiplying Eq. (|13(l by each of the 
functions bi{mj.), one obtains K equations. The integra- 
tion over the whole volume of the star, V, of each of these 
K equations gives 

/ pr^^b.dV = [ div {pr^D'grsLd{n))bidV. (15) 

Jv Jy 

Using the general relations: 

div(av) — adiv(v) + grad(a) • v, (16) 

div(v)dF= / vdS, (17) 
V Js 

where a is a scalar, v a vector and S the surface of the 
volume V, one obtains that the right hand term of Eq. (|15|l 
becomes 



/ pr^ D' grad{n)b, dS- j pr^ D' grad{n) ■ grad{b^) dV. 
Js Jv 

The integral on the surface is null, since grad(ri) = on 
the surface. Thus one has K equations of the type 



pr^^b,dV+ I pr'D'^^dV = 0. 

or or 



(18) 



Using drrir ~ pdV = Anr^pdr, one obtains 



M 



^D'{^Tir^pf 



on db, 
drrir dnir 



drrir = 0. (19) 



The functions bi constitute a set of independent functions, 
therefore the unknown function 51 (to^, t) cai\ be expressed 
as a linear combination of 6.;'s. One can write 



K 



Q,(mr, t) — fljbj(mr), 



(20) 



where ilj — il{mj,t). Equation H20|l simply says that to 
obtain fl{mr,t), with rrir G [wi+i,mi], one simply inter- 
polates linearly as a function of mass between f2i+i and 

Expressing Q{mr,t) as indicated in Eq. (PUjl . Eq. l(Tn|l 

becomes 



J 

where 



dt 



Aij — / r^bibj drrir 
Jo 



AI 



B,. = / r^D'{A^r'pf^^ dm. 

Orrir Orrir 



(21) 



(22) 



(23) 



One has used the facts that the fij's do not depend on 
TUr and the &j's do not depend on time. When |j — i| > 1, 
there is no overlap between the functions bi and bj, thus 



Aij = 0, for \j -i\>l. 



(24) 



different from zero. The same is true for the Bij. Let us 



For a given value of i only, Ai^i-i, Ai^i and Ai^i^i are 
different fro 
estimate Ai 

Ai^i-i 



i.i— 1 • 



rrir — rrii Jrir 



mi 



mi_i - rrii rrii - "^i-i 



dm., 



(25) 



one approximates by 0.5 (r^(mi_i) + r^(mi)), then, by 
trivial integration, one obtains 



A, 



-0.5 (r2(m,_i) -f r2(m,)) 



rrii ~ rrii-i 
6 



(26) 



The other matrix elements are obtained in the same way: 

(27) 

A,,, = 2(A,,,_i+A:+i,0. (28) 



A+i, = -0.5 (r2(m,+i) + r\m,)) !!h±l_^ 



At the centre and at the surface, one has 



Ak,k 
Ak,k-i 

-42,1 

^1,1 



Q.->r\mK-i)^, 

O.SAk.k, 

0.5Ai,i, 

0.5 (r2(l)-Hr2(2)) 
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Let us set D — {iirr^ p)^ D' and Ui^i-i an appropriate 
mean value of D between the shells i and i—1 (see Sect . 4) , 

(29) 



then the Bij becomes 



B 



D 



B 



nii ' 



(30) 
(31) 



Bi^i = —Bi^i-i — Bi^i^i. 
At the centre and at the surface, one obtains 

Bk,k-i = —Bk,k, 
^2,1 = —Bi^i, 

^■^ mi— m2 

Adopting an implicit discretisation in time, we have for 
Eq. ^ 



71+1 



At 



Y,B,,n]+' = o, 



or 



E 



At 



B. 



E 



A. 



At ^ 



(32) 



(33) 



where £7" is equal to il{mj, i"). This linear system of equa- 
tions is then solved by standard procedure. 

3.2. The Crank-Nicholson finite differences method 

Let us consider in a first step the case of the diffusion of 
the chemical elements. Multiplying Eq. Hl()|l by Anr"^ and 
integrating with respect to the spatial coordinate r from 
r — Ar/2 to r + Ar/2, one obtains 



Anpr D- 



where Am — ^T^r^pdr. The quantity Ar has been 



dr 



r+Ar/2 
r-Ar/2 



(34) 



-Ar/2 



chosen sufficiently small for to remain constant over 
the spatial range Ar around r. This integration reduces 
the equation to a first order differential equation, which 
simply expresses the fact that the time variation of the 
mass of element ? in a spherical shell is equal to the dif- 
ference between the mass of element i which enters in the 
shell and the mass of this same element which goes out 
from the shell. As indicated in Sect. 2, one can interpret 
the quantity, 




r,=0 r,_ 



'kti 'k 'k-i 
Stellar radius 



Fig. 1. Discretized distribution of the abundance of ele- 
ment I in a star model as a function of radius. Only a 
few mesh points are represented. K is the total number of 
mesh points. Numbering increases from outside to inside. 

as a diffusion velocity, < >, for the element i (let us 
recall that the minus sign appears because, when the spa- 
tial abundance gradient is positive, the velocity is directed 
inward). In that case, Eq. H34|l writes 



Am 



dt 



-Anpr^Xi <Vi> 



|r-|-Ar/2 
lr-Ar/2' 



which makes clearly appear the difference of the mass 
fluxes between the two borders of the shell. This is a more 
elementary description of diffusion. 

One can also derive Eq. H35|l from the continuity equa- 
tion. Indeed supposing that the bulk gas velocity is zero, 
that dp/dt = (stationary situation), and that there is 
no sink/source terms of matter, the continuity equation 
becomes 



P- 



dX^ 
dt 



j,2 Qj, 



{pr^ <V,>X, 



(36) 



By equating the right hand sides of Eqs. H10() and (|36() . 
one obtains for < Vi > Eq. H35|) above. 

Let us now consider the discretized abundance profile 
of element i sketched in Fig. 1. The mass of element i 
removed from point fc + 1 and added to point k per unit 
time is (see Eg. I34|l : 



"k+l/2Pk+l/2Dk+l/2 



dX, 
dr 



(37) 



fe+l/2 



(35) 



where the physical quantities with a lower subscript fc-|-l/2 
are estimated at the midpoint between the mesh points 
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k+1 and k. For instance, one takes 7'fc+i/2 — 0-5{rk+rk+i)- 
For the diflFusion coefficient, we use tlie procedure exposed 
in Sect. 4. A similar expression can be written for the mass 
of element i removed from point k and added to point k—1. 
Let us set Xi ^, the mass fraction of element i evaluated 
at the mesh point k. Thus one has. 



Arrifc 



dt 



A 2 dXi 
-47rr;,+i/2Pfc+i/2^fe+i/2-^ 



-47rr^_;^/2P/c-i/2£'fc~i/2 



dX, 
dr 



k+l/2 



fe-1/2 



(38) 



where Amt = (^'^ ^''^ A.-kt^ pdr. Let us now discretize this 

•"'fc + l/2 ' 

equation as a function of time. The left hand term of 
Eq. H38|l can be written 



Awfe 



d{X,,k) 
dt 



Aruk- 



X. 



n+l 
i.k 



X] 



(39) 



where the superscript n indicates that the quantity is esti- 
mated at the time t" . The Crank-Nicholson method con- 
sists in evaluating the right hand term at the same time 
as the left hand term, i.e. at time t"+^/^. The system is 
centered in time and thus second order accurate in time. 
One obtains 



At 

O'fc+l/2 

Amfc 

ffc-l/2 



Dk+1/2 — 



ffe+1 - ffe 



^i^k^+^i'.k ^?,k-l+^zk 



Dk-1/2 — 



. ^ , (40) 

where one has introduced the quantity o-k^xji defined by 

c^/c+i/2 = 47rr^^;^/2Pfc+i/2- (41) 
Let us define the quantity [A; + 1 , fc] by 



1 gfc+l/2At 

\k -f 1, A:J = -0-/C+1/2 , 

2 Tfc+i - rfc 



(42) 



and let us separate the abundances evaluated at time 
from those estimated at time Eq. (|40|l then becomes 



[fc, k 1] „+i 



Avuh 



^"^1 + 1 



[/c,fc-l] [fc+l,fc] 
Amfc ArTifc 



X 



n+l 

i,k 



[k + 1, k] 



[k,k-l] 

Anik 



X 



- ^ A ^ [k,k-l] [k + l,k] \ ^ 

+ ^ A ' A X 

\ Amfc Amfc / 

[k + l,k 



Amfc ^''=+1 



(43) 



Following a similar line of reasoning, one obtains for the 
equations at the center and at the surface (see Fig.^, 



ArriK 
ArriK 



X 



i.K-l 



[A', K-l] 
AniK 

ArriK 



XI 



K 



X 



i.K 



(44) 



[2,1] _ 
Ami 



X 



i,2 



Ami 



1 - 



[2,1] 
Ami 

[2,1 



Ami 



vn+l _ 
^i.l - 



(45) 



where ArtiK — Jq"^ ^''^ Anr^pdr and Ami = 
■I'ri 1/2 ^^''"^ Pd'"'" ■ Eqs. 63, ESIl and constitute a sys- 
tem of linear equations whose unknowns are the X"'^^ . It 
is solved by using classical methods of tridiagonal matrix 
inversion. 

The discretized equations describing the diffusion of 
the angular momentum are obtained in a similar way. The 
final result are equations similar to Eqs. (|43|1 (|44|l and 
(I45|l with Xi replaced by fi, Amfc replaced by Amfcr^., 
crfe+i/2 replaced by o-fc+i/2'^fc+i/2 and CTfc_i/2 replaced by 
o'fe-i/2^fc_i/2- Of course the diffusion coefficient must be 
the one describing the diffusion of the angular momentum. 

Let us emphasize here that this method is not fully im- 
plicit and thus, one can expect that some sort of Courant's 
condition will limit its domain of validity (see Sect. 5). 

3.3. The implicit finite differences metiiod 

In the implicit finite differences method, the right hand 
term of Eq. (|38|l is estimated at time t""*"^ (in the explicit 
method the right hand term would be estimated at time 
i"). In this case, one obtains 



^i.k - ^Lk 



At 

O'fe+l/2 

Amu 

ffe-l/2 



Dk+1/: 



Dk-i/: 



rk+1 - Tk 

^i.k ^i.k-l 



(46) 



Amfc ' Tfc - rfc_i 

Let us define the bracket terms [k + l,k] by 

Jfc.+ l/2At 
fc + l,fc = crfc+i/2 , 

Tk+l - Tk 

doing so, we obtain, separating the abundances at time t" 
from those estimated at time 



[fc,fc-l] ,.n+l [fe, fc - 1] _ [k±lM\ n+l 

Amfc '^'^-i I Amfc Amfc j 



[fc -f 1, fc] ^n+l _ 

Amfc ^^.'=+1 - ^^.'^ 



(47) 
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Following a similar line of reasoning, one obtains for the 
equations at the center and at the surface, 



[K,K-1] 



vn+l 



[2,1] 
Ami 



X 



n+l 



1 - 



1 - 



[2,1] 



[K,K-1] 



= Kk , (48) 



A' 



mi 



X 



n+l 



= x. 



i.l 



(49) 



One can easily check that the system of discrctizcd 
equations (both in the cases of the Crank-Nicholson 
scheme and in the implicit finite differences method) con- 
serves the integrated quantities of the elements over the 
mass of the star. The system also keeps equal to one the 
sum of the X^ and does not induce any diffusion when the 
chemical gradients are flat. 

For obtaining the equations for the angular momen- 
tum, we apply the same recipe as indicated at the end of 
the previous section. The system of equations describing 
the transport of the angular momentum also conserves the 
total angular momentum. 

4. Estimate of 

Let us consider two mesh points as in Fig. 13 Let us sup- 
pose that the diffusion coefficient D^-i is the coefficient 
for the whole region between rfc_i and r^. Similarly the 
diffusion coefficient Dk is the coefficient for the whole re- 
gion between Tc and r^. When both mesh points are in 
a radiative or convective region, the radius Tc is equal to 
{fk-i + "fk)l'^i otherwise r^. is taken as the radius of the 
limit between the radiative and the convective zone. 

Over the path /Ar (see Fig.lJJ, the particles have an 
average diffusive velocity V/c (see Ea. lH5|l and over the path 
(1 — /)Ar they have a diffusive velocity < T4_i >. The 
total time for going from fc to /c — 1 is equal to 



At = /Ar/ < 14 > +(1 - /)Ar/ < 14,_i > . 



(50) 



The mean velocity < Vfc_i/2 > over the whole interval is 
given by Ar/At, thus one has 



< yfc_i/2 >= 



1 



j_ 



(1-/) 

<Vfc> ' <yfe_i> 
^ < Vk >< Vk-i > 

f < Vk-i > +(1 - /) < > ■ 

Now, the diffusive coefficient is proportional to the diffu- 
sive velocity (see Eg. This suggests the way to com- 
pute the diffusive coefficient between two shells: 



D 



k~l/2 



Dk-iDk 



fDk-i + {l- f)Dk 
This expression is more physical than 



(51) 



fe-l/2 



fe-i 



AO/2, 

that simply averages the diffusion coefficients. 
Equation l(5T|l implies that if Dk » Dk-i and (1 — /) is 



fAr 



D,. 



(l-f)Ar 



D„ 



- k-l 



Stellar radius 

Fig. 2. Schematic representation of the mass shell between 
the radii and rk-i- The diffusive coefficient Dk oper- 
ates in the zone between and Tc- The diffusive coeffi- 
cient Dk-i operates in the zone between Tc and rk~i. The 
quantity Ar is equal to rfc_i — rk- 



of the order of 0.5 then Dk-1/2 ^ Dk-i- As expected, the 
smallest diffusion coefficient governs the diffusion between 
the two mesh points. The simple algebraic mean would 
give Dk-1/2 ^ Dk/2 a much greater diffusion coefficient, 
which is not physically justified. One sees also that when 



Dk 



fe-i 



Dk 



D then D 



fc-1/2 



D whatever the value 



of /. In Sect. 5 below we shall illustrate by a numerical 
example the importance of correctly evaluating the 
diffusion coefficient at the interface between a convective 
and a radiative zone. 

When the finite elements method is used, we need to 
evaluate the mean value between two mass shells of the 
diffusion coefficient, D, multiplied by a factor a, which 
for the transport of the chemical elements is (47rr^p)^ and 
for that of the angular momentum is r^(47rr^p)^ (see Sect. 
3.1). In that case, one adopts the following procedure: 



afc-1/2 A-1/2 = 0.5 (ttfc -I- ttfc+i) 



Dk^iDk 



fDk 



{l-f)Dk 



5. Tests and comparisons 

5.1. Initial configuration and Cou rant's condition 

In this section, we study how the different methods de- 
scribed above behave in a very simple configuration. Let 
us consider a uniform density sphere of 1 Mq, with a ra- 
dius of 1 R0, composed of two chemical elements X and 
Y. The initial distribution of these two elements in the 
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sphere is a quasi step function. The variation as a func- 
tion of the radius of the abundance of one of these two 
elements, let us call it X, at the beginning of the compu- 
tation is shown in Fig. 13 by the dotted line, (cf. also the 
dotted hues in Figs.|Hlto|SJ|. At the center, its initial value 
is 0.99, in the outer regions X = 0.01. The abundance Y is 
simply 1 — X. The diffusion coefficient is set equal to 10^* 
cni^ s^^ in the interior region {i.e. for a radius r inferior to 
0.7963 R© or a mass inferior to 0.5049 Mq) and to 1.5 10^ 
cm^ s~^ in the outer zone. Such values for the diffusion 
coefficient imply that the interior has a very small mixing 
timescale (of the order of B? /D ~ 1 yr) and therefore will 
be always homogenized by mixing, while the outer one 
will have a much longer mixing timescale (of the order of 
4.26 Myr). Such a situation is quite analog to a convective 
core surrounded by a radiative envelope in a star. In the 
following we shall define the interior region as "the core" , 
and the outer one as "the envelope" . In our numerical ex- 
periments we keep the diffusion coefficient constant as a 
function of time. The sphere is decomposed in 100 shells 
of equal mass. 

From these data one expects that the whole star will be 
completely mixed in less than about 10^ yr and that the 
final abundance of element X in the homogeneous star will 
be 0.5048. From this initial structure, on can also estimate 
the "Courant condition" . Let us recall that the "Courant 
condition" imposes a superior limit Aic to the time step 
for an explicit method to be stable (see e.g. Press et al. 
I1992|l . To estimate this limit, one has to compute for each 
element and for each mass shell, the time required for the 
element to diffuse through the width of the shell i.e. 



Ar 



Ar 



VY R\^\ 
X \ dr \ 



Then one has to take the smallest value of all these times. 
In a discretized form, Ate may be written 



Atr = Min 



in 



D.. 



+ 1/2 



Xi-X; 



+ 1 



For the simple initial structure considered here, the 
diffusion time through a shell takes a non infinite value 
only at the interface between the core and the envelope. 
In this case, one has that — r^+i is equal to 0.000199 
Rq, £'i+i/2j estimated from Eq. (|51l) with / taken equal 
to 0.5, is equal to 3 10^ cm^ s, X,+i - X^ ^ 0.98 and 
x.-i-x.+i gqual to 0.5. This gives 

Ate ^ 1 yr. 

Let us stress that implicit methods are generally stable 
for any time step, and thus arc not limited by the Courant 
condition exposed above (cf. Press et al. I1992|l . In the 
following we shall test this point. 



Table 1. Maximum value over the star of the quantity 
Ax = 1 - X[k) - Y{k) and of AB = ^""1^7.^"""' , for 
different values of r. At (in years) and for different nu- 
merical schemes. X and Y are the mass fraction of the 
two elements composing the "star", B is the total angu- 
lar momentum of the star. The labels CN and FI are for 
"Crank-Nicholson" and "Fully Implicit" respectively (see 
text). 
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-484 
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29542 


2.24 


3772 


-0.67 


12426 


-0.80 
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245311 


-1.34 


37871 


0.44 


123274 


-0.90 


10" 


10^ 


245311 


-1.72 


25155 


-0.65 


48542 


-0.80 



5.2. Comparisons of the methods 

Let us begin by illustrating the importance of correctly es- 
timating the diffusion coefficient at the interface between 
a convective and a radiative zone. On Fig. 13 the result mg 
distributions of the element X inside the star after 1000 yr 
is shown. The time step used is At ~ 10 yr The continuous 
line shows the results obtained using Eq. (|51|) for £'fc_i/2- 
The dashed line represent the solution obtained using a 
simple algebraic mean. We see that the results are signif- 
icantly different, in particular the algebraic mean tends 
to slightly increase the convective core. This may have 
important consequences when such increases are repeat- 
edly applied over the whole evolution of a star. The use of 
Eq. (|51|l which results from physical considerations is thus 
recommended, and this the one we used in the numerical 
experiments we now describe. 

We have computed the diffusion of the chemical ele- 
ments (and of the angular momentum, see Sect. 5.7) for 
different durations t of the period during which the dif- 
fusion operates and for different time steps At. In Fig. ^ 
the set of values (log At, logr) explored are indicated. For 
each couple of values, we performed computations with the 
three methods described above, namely the implicit finite 
elements method, the Crank-Nicholson finite differences 
method and the implicit finite differences method. Most 
of the results are displayed in Figs. |5l to |S| On these fig- 
ures, the results obtained for increasing durations, using 
the same time step, are ordered horizontally (from left to 
right), while the results obtained for the same duration, 
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0.795 0.8 0.805 

Fig. 3. Comparison of the abundances obtained with the 
imphcit finite differences method using two different ways 
of evaluating the diffusion coefficient -Dfc-1/2 between two 
mesh points (see Sect. 4). The dotted line shows the initial 
situation [t—Q). A time step of 10 years was adopted 
and the computation was stopped at an age of r = 1000 
years. The dashed line show the result obtained when one 
takes for -Dfc_i/2 a simple algebraic mean, the continuous 
line presents the result obtained using Eq. H51|l . 



with increasing time steps, are arranged vertically (from 
top to bottom). Since the Courant time step is equal to 
one year, the time step At expressed in years gives directly 
the time step in units of the Courant time step. 
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Fig. 4. Computations have been performed with the three 
methods for each couple of (log At, logr) values corre- 
sponding to the position of a square in the above figure; 
At is the time step, t is the duration over which the com- 
putation was performed, (no computation, of course, have 
been performed in the grey zone where log At > log t) . 
The zone I (square with a cross inside) shows the "non- 
physical" region for the implicit finite elements method. 
By "non-physical" , we mean here that negative abun- 
dances are obtained. The zone II (square with a filled tri- 
angle inside) corresponds to the zone of "non-physical" so- 
lutions for the Crank-Nicholson finite differences method. 
The implicit finite differences method gives physical solu- 
tion in all the cases considered here (indicated by squares 
filled or empty). 



5.3. The implicit finite elements method 

The solutions obtained with the implicit finite elements 
method are shown in Figs. |31 and El (the dashed-dotted 
lines) . This method gives reasonable results in all the cases 
explored in this work except for small values of r. We 
note that for r < 1000 yr, whatever is the time step, 
the solution leads to negative values of X just above the 
"convective core" . If the integration is performed over a 
sufficiently long time, the system eventually reaches a rea- 
sonable solution, even when instabilities appear at earlier 
time (see for instance in Figs. and the evolution when 
r increases with At = 0.5 yr). One notes that the star 
is completely mixed for r > lO'', i.e. for durations more 
than twice the mixing timescale of the envelope (about 
4.3 Myr, see Sect. 5.1), which is not very satisfactory. The 
final abundance in the homogeneous star is on the other 
hand equal to that expected (~ 0.5). 

As expected from an implicit scheme, reasonable so- 
lutions are obtained even if the time steps are much 



greater than the time step given by the Courant condition. 
Inspection of the results for r > 1000 yr show in general a 
great stability of the solution with respect to the choice of 
the time step. The solutions obtained with the implicit fi- 
nite elements method are in general less mixed that those 
obtained with the implicit finite differences method (see 
Figs. □ and EJ. 

5.4. The Crank-Nicholson finite differences method 

This method is a kind of compromise between a fully im- 
plicit and an explicit scheme. In order to better under- 
stand what happens here, let us briefly recall a few gen- 
eralities about implicit and explicit schemes (see Press et 
alflW, p. 838-842). 

Explicit finite differences schemes are stable only if the 
time steps satisfy the Courant condition. Such methods 
are not suited for the computation of the secular evolu- 
tion of stars. Indeed, we are interested in modeling the 
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evolution of features with spatial scales of the order of the 
radius of the star R, much greater than the distance be- 
tween two mesh points Ar. If we are limited to time steps 
satisfying the Courant condition, we will need of the or- 
der of i?^/(Ar)^ steps before anything interesting begins 
to happen. We would like instead use timesteps of the or- 
der of B? /D or, maybe, for purpose of accuracy, somewhat 
smaller. 

With such great timesteps, it is no long possible to 
describe accurately what happens at small spatial scales. 
However at small scale, the differencing scheme must do 
"something stable, innocuous, and perhaps not too phys- 
ically unreasonable" write Press et al. H1992fl . One pos- 
sibility is to use a Crank-Nicholson differencing scheme. 
One of its main property is to let small-spatial-scale fea- 
tures maintain their initial amplitudes. In that 
cording to Press et al. H1992|l . the evolution of the larger- 
scale features, in which we are interested in, take place su- 
perposed with a kind of "frozen in" (though fluctuating) 
background of small-scale features. This is what happens 
in our numerical experiments. Indeed looking at the con- 
tinuous Hnes in Figs. El and one sees that the amplitude 
of the initial step in chemical composition is more or less 
maintained, although with great fluctuations. 

The method gives reasonable solutions for not too big 
time steps. More precisely, physical solutions are obtained 
when Ai in years is inferior to lO-^/r, or when r in years 
is inferior to 100 yr (see Fig. 2)). This restriction on the 
time step is reminiscent of a kind of Courant 's condition. 
This is not so surprising given that the Crank-Nicholson 
scheme is a mixture of both an implicit (not submitted to 
Courant condition) and an explicit method (submitted to 
Courant condition). 

Let us finally note that in its domain of validity in the 
log At versus log r plane this method gives the same solu- 
tion as the implicit finite differences method, and since, it 
is centered in time, this scheme is second-order accurate 
in time. 

5.5. The implicit finite differences method 

Another possibility for imposing an "inoccuous" be- 
haviour to the small-scale features is to adopt an implicit 
method. Such schemes drive small-scale features to their 
equilibrium form, i.e. imposes that at small scales 



d_ 

dr 



pr^D- 



dr 



0. 



This can be seen from Eq. 146|l with At oo. 

The solutions obtained with the implicit finite differ- 
ences method are shown in Figs. 13 and |H1 (continuous 
lines). The results show in general a great stability of the 
solution with respect to the choice of the time step. Only 
when the time step is of the same order of magnitude as r 
can we notice some differences (compare for instance the 
continuous line for r = 10^ yr and At — 10^ yr with the 
case T = 10"^ yr and At = 10'' yr in Fig.lHJ. 



This method leads to a completely mixed star after a 
time less than 10^ in agreement with the estimate made 
in Sect. 5.1. The final abundance of the element X in the 
homogeneous star is also equal to that expected. 

Although, in that case, the differencing scheme is only 
first-order accurate in time, it has the advantage over the 
Crank-Nicholson finite differences method and the im- 
plicit finite elements method to propose a physical solu- 
tion {i.e. without negative abundance values) for all the 
cases investigated here (see Fig.QJ. Moreover, as seen just 
above, it predicts a timescale for the star to be completely 
homogenized in agreement with the usual analytical esti- 
mate (see Sect. 5.1). In that respect this method does 
appear as the best one. 

5.6. Conservation of the sum of the abundances 

Let us now have a look on the ability of the different 
schemes to conserve the sum of the abundances (in mass 
fraction). In each shell and at any time step, one should 
have that X+Y^l. In Tabled we indicate for all the cases 
where the three methods give physical solutions, the max- 
imum over the star of the quantity Ax = 1 — X{k) — Y(k) 
at the end of the computation. As expected we see that in 
general A^ becomes greater when r increases. The Crank- 
Nicholson finite differences method appears to give, in 
most of the cases, the smallest values for Ax (inferior 
to 10"'' for the cases considered). This is likely related 
to the fact that this method is second order accurate in 
time. The implicit finite elements method gives in general 
the highest values (in the worst case Ax is of the order of 
10"'^), while the implicit finite differences method gives in 
general results inbetween. Thus the implicit finite differ- 
ences method enables to avoid unphysical solutions and 
keeps reasonably well the sum of the abundances equal to 
one. 

5.7. Conservation of the angular momentum 

We have performed similar tests and comparisons for the 
diffusion of the angular momentum. We started from an 
initial configuration where the "convective core" defined 
in Sect 5.1 has an angular velocity equal to 1-10~^ sec~^ 
and the radiative envelope an angular velocity equal to 
1T0~^ sec~^. The other variables were taken as described 
in Sect. 5.1. In such situation the Courant time step, which 
writes 



At^ = Min 



in 



n+iY 



D 



i+l/2 



is equal to 1.24 year, not very different from the one we 
obtained for the diffusion of the chemical species. It is 
thus not surprising that the results we obtained are qual- 
itatively similar to those presented in Figs. El to |H1 The 
zones where unphysical solution are encountered are the 
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same as those presented in Fig.0] We show in Tabled the 
value of the quantity AB which is equal to 

A j-f -^final -^initial 

Ai( = , 

-L^initial 

where i?initiai, ^finai are the total angular momentum of 
the star at the beginning, respectively at the end of the 
computation. The relative error on the total angular mo- 
mentum is much greater than the error on the sum of the 
abundances. This arises because the total angular momen- 
tum is an integrated value over the whole system, imply- 
ing that the errors can accumulate. The sum of the abun- 
dances, instead, is evaluated locally in a given shell. 

We see that the finite differences method appear to 
better conserve the total angular momentum. The Crank- 
Nicholson finite differences method gives the best results, 
followed by the implicit finite differences method which 
reaches the same level of accuracy in most cases. Thus the 
implicit finite differences method that we recommended 
on the base of the results obtained for the diffusion of the 
chemical species gives also very satisfactory results for the 
diffusion of the angular momentum. 
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6. Conclusion 

From the above numerical experiments, it appears that 
the implicit finite differences method seems to be the most 
robust one, giving physical solutions in all the cases stud- 
ied here spanning more than nine orders of magnitude in 
T and At. It has moreover the following characteristics: 
1) it reduces the problem to a first order differential equa- 
tion, 2) it enables an easy and clear interpretation of what 
happens physically in the system, 3) it is quite easy to im- 
plement in a code, 4) it conserves reasonably well the sum 
of the abundances at each mesh point as well as the total 
angular momentum. 

On the basis of the new tests and analysis made here, 
we can recommend this method to resolve the diffusion 
equation in stellar interiors. Of course many more tests 
could be performed by changing the initial conditions. We 
restrained our discussion here on a case with a very sharp 
gradient in order to test the different numerical methods in 
some extreme conditions. Adopting an initially shallower 
gradients will tend to make the things much more easier 
for all three methods and they would give identical results. 

The diffusion coefficient between two mesh points has 
to be evaluated correctly in order to obtain reliable re- 
sults. The mean diffusion coefficient is equal neither to 
the algebraic nor to the geometric mean. Its expression is 
given in Eq. ()51|l . 
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Fig. 5. Abundance of the element X as a function of the distance to the center. The dotted line show the situation at 
the beginning of the computation. The continuous line refers to the solution obtained after a time r by the Crank- 
Nicholson finite differences method. The dashed-dotted line shows the result obtained after a time r by the implicit 
finite elements method. The time step At used in each case is indicated. In all the cases, the Courant condition time 
is about 1 yr. This means also that the value of At is equal to the ratio At/ Ate- 
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Fig. 6. Same as Fig. 5 for other values of r and At. 
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Fig. 7. Same as Fig. 5 except that the continuous line corresponds to the result obtained after a time r by the implicit 
finite differences method. 
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Fig. 8. Same as Fig. 6 except that the continuous line corresponds to the result obtained after a time r by the implicit 
finite differences method. 



